C
C     ..................................................................
C
C        SUBROUTINE DPRQD
C
C        PURPOSE
C           CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN POLYNOMIAL
C           WITH REAL COEFFICIENTS.
C
C        USAGE
C           CALL DPRQD(C,IC,Q,E,POL,IR,IER)
C
C        DESCRIPTION OF PARAMETERS
C           C     - COEFFICIENT VECTOR OF GIVEN POLYNOMIAL
C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
C                   THE GIVEN COEFFICIENT VECTOR GETS DIVIDED BY THE
C                   LAST NONZERO TERM
C                   DOUBLE PRECISION ARRAY
C           IC    - DIMENSION OF VECTOR C
C           Q     - WORKING STORAGE OF DIMENSION IC
C                   ON RETURN Q CONTAINS REAL PARTS OF ROOTS
C                   DOUBLE PRECISION ARRAY
C           E     - WORKING STORAGE OF DIMENSION IC
C                   ON RETURN E CONTAINS COMPLEX PARTS OF ROOTS
C                   DOUBLE PRECISION ARRAY
C           POL   - WORKING STORAGE OF DIMENSION IC
C                   ON RETURN POL CONTAINS THE COEFFICIENTS OF THE
C                   POLYNOMIAL WITH CALCULATED ROOTS
C                   THIS RESULTING COEFFICIENT VECTOR HAS DIMENSION IR+1
C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
C                   DOUBLE PRECISION ARRAY
C           IR    - NUMBER OF CALCULATED ROOTS
C                   NORMALLY IR IS EQUAL TO DIMENSION IC MINUS ONE
C           IER   - RESULTING ERROR PARAMETER. SEE REMARKS
C
C        REMARKS
C           THE REAL PART OF THE ROOTS IS STORED IN Q(1) UP TO Q(IR)
C           CORRESPONDING COMPLEX PARTS ARE STORED IN E(1) UP TO E(IR).
C           IER = 0 MEANS NO ERRORS
C           IER = 1 MEANS NO CONVERGENCE WITH FEASIBLE TOLERANCE
C           IER = 2 MEANS POLYNOMIAL IS DEGENERATE (CONSTANT OR ZERO)
C           IER = 3 MEANS SUBROUTINE WAS ABANDONED DUE TO ZERO DIVISOR
C           IER = 4 MEANS THERE EXISTS NO S-FRACTION
C           IER =-1 MEANS CALCULATED COEFFICIENT VECTOR REVEALS POOR
C                   ACCURACY OF THE CALCULATED ROOTS.
C                   THE CALCULATED COEFFICIENT VECTOR HAS LESS THAN
C                   6 CORRECT DIGITS.
C           THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED
C           COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE BEEN
C           CALCULATED.
C           THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT VECTOR IS
C           RECORDED IN Q(IR+1).
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF
C           THE QUOTIENT-DIFFERENCE ALGORITHM WITH DISPLACEMENT.
C           REFERENCE
C           H.RUTISHAUSER, DER QUOTIENTEN-DIFFERENZEN-ALGORITHMUS,
C           BIRKHAEUSER, BASEL/STUTTGART, 1957.
C
C     ..................................................................
C
      SUBROUTINE DPRQD(C,IC,Q,E,POL,IR,IER)
C
C      DIMENSIONED DUMMY VARIABLES
      DIMENSION E(1),Q(1),C(1),POL(1)
      DOUBLE PRECISION Q,E,O,P,T,EXPT,ESAV,U,V,W,C,POL,EPS
C
C        NORMALIZATION OF GIVEN POLYNOMIAL
C           TEST OF DIMENSION
C        IR CONTAINS INDEX OF HIGHEST COEFFICIENT
      IR=IC
      IER=0
      EPS=1.D-16
      TOL=1.E-6
      LIMIT=10*IC
      KOUNT=0
    1 IF(IR-1)79,79,2
C
C        DROP TRAILING ZERO COEFFICIENTS
    2 IF(C(IR))4,3,4
    3 IR=IR-1
      GOTO 1
C
C           REARRANGEMENT OF GIVEN POLYNOMIAL
C        EXTRACTION OF ZERO ROOTS
    4 O=1.0D0/C(IR)
      IEND=IR-1
      ISTA=1
      NSAV=IR+1
      JBEG=1
C
C        Q(J)=1.
C        Q(J+I)=C(IR-I)/C(IR)
C        Q(IR)=C(J)/C(IR)
C        WHERE J IS THE INDEX OF THE LOWEST NONZERO COEFFICIENT
      DO 9 I=1,IR
      J=NSAV-I
      IF(C(I))7,5,7
    5 GOTO(6,8),JBEG
    6 NSAV=NSAV+1
      Q(ISTA)=0.D0
      E(ISTA)=0.D0
      ISTA=ISTA+1
      GOTO 9
    7 JBEG=2
    8 Q(J)=C(I)*O
      C(I)=Q(J)
    9 CONTINUE
C
C           INITIALIZATION
      ESAV=0.D0
      Q(ISTA)=0.D0
   10 NSAV=IR
C
C        COMPUTATION OF DERIVATIVE
      EXPT=IR-ISTA
      E(ISTA)=EXPT
      DO 11 I=ISTA,IEND
      EXPT=EXPT-1.0D0
      POL(I+1)=EPS*DABS(Q(I+1))+EPS
   11 E(I+1)=Q(I+1)*EXPT
C
C        TEST OF REMAINING DIMENSION
      IF(ISTA-IEND)12,20,60
   12 JEND=IEND-1
C
C        COMPUTATION OF S-FRACTION
      DO 19 I=ISTA,JEND
      IF(I-ISTA)13,16,13
   13 IF(DABS(E(I))-POL(I+1))14,14,16
C
C        THE GIVEN POLYNOMIAL HAS MULTIPLE ROOTS, THE COEFFICIENTS OF
C        THE COMMON FACTOR ARE STORED FROM Q(NSAV) UP TO Q(IR)
   14 NSAV=I
      DO 15 K=I,JEND
      IF(DABS(E(K))-POL(K+1))15,15,80
   15 CONTINUE
      GOTO 21
C
C           EUCLIDEAN ALGORITHM
   16 DO 19 K=I,IEND
      E(K+1)=E(K+1)/E(I)
      Q(K+1)=E(K+1)-Q(K+1)
      IF(K-I)18,17,18
C
C        TEST FOR SMALL DIVISOR
   17 IF(DABS(Q(I+1))-POL(I+1))80,80,19
   18 Q(K+1)=Q(K+1)/Q(I+1)
      POL(K+1)=POL(K+1)/DABS(Q(I+1))
      E(K)=Q(K+1)-E(K)
   19 CONTINUE
   20 Q(IR)=-Q(IR)
C
C           THE DISPLACEMENT EXPT IS SET TO 0 AUTOMATICALLY.
C           E(ISTA)=0.,Q(ISTA+1),...,E(NSAV-1),Q(NSAV),E(NSAV)=0.,
C           FORM A DIAGONAL OF THE QD-ARRAY.
C        INITIALIZATION OF BOUNDARY VALUES
   21 E(ISTA)=0.D0
      NRAN=NSAV-1
   22 E(NRAN+1)=0.D0
C
C           TEST FOR LINEAR OR CONSTANT FACTOR
C        NRAN-ISTA IS DEGREE-1
      IF(NRAN-ISTA)24,23,31
C
C        LINEAR FACTOR
   23 Q(ISTA+1)=Q(ISTA+1)+EXPT
      E(ISTA+1)=0.D0
C
C        TEST FOR UNFACTORED COMMON DIVISOR
   24 E(ISTA)=ESAV
      IF(IR-NSAV)60,60,25
C
C        INITIALIZE QD-ALGORITHM FOR COMMON DIVISOR
   25 ISTA=NSAV
      ESAV=E(ISTA)
      GOTO 10
C
C        COMPUTATION OF ROOT PAIR
   26 P=P+EXPT
C
C        TEST FOR REALITY
      IF(O)27,28,28
C
C        COMPLEX ROOT PAIR
   27 Q(NRAN)=P
      Q(NRAN+1)=P
      E(NRAN)=T
      E(NRAN+1)=-T
      GOTO 29
C
C        REAL ROOT PAIR
   28 Q(NRAN)=P-T
      Q(NRAN+1)=P+T
      E(NRAN)=0.D0
C
C           REDUCTION OF DEGREE BY 2 (DEFLATION)
   29 NRAN=NRAN-2
      GOTO 22
C
C        COMPUTATION OF REAL ROOT
   30 Q(NRAN+1)=EXPT+P
C
C           REDUCTION OF DEGREE BY 1 (DEFLATION)
      NRAN=NRAN-1
      GOTO 22
C
C        START QD-ITERATION
   31 JBEG=ISTA+1
      JEND=NRAN-1
      TEPS=EPS
      TDELT=1.E-2
   32 KOUNT=KOUNT+1
      P=Q(NRAN+1)
      R=ABS(SNGL(E(NRAN)))
C
C           TEST FOR CONVERGENCE
      IF(R-TEPS)30,30,33
   33 S=ABS(SNGL(E(JEND)))
C
C        IS THERE A REAL ROOT NEXT
      IF(S-R)38,38,34
C
C        IS DISPLACEMENT SMALL ENOUGH
   34 IF(R-TDELT)36,35,35
   35 P=0.D0
   36 O=P
      DO 37 J=JBEG,NRAN
      Q(J)=Q(J)+E(J)-E(J-1)-O
C
C           TEST FOR SMALL DIVISOR
      IF(DABS(Q(J))-POL(J))81,81,37
   37 E(J)=Q(J+1)*E(J)/Q(J)
      Q(NRAN+1)=-E(NRAN)+Q(NRAN+1)-O
      GOTO 54
C
C        CALCULATE DISPLACEMENT FOR DOUBLE ROOTS
C           QUADRATIC EQUATION FOR DOUBLE ROOTS
C           X**2-(Q(NRAN)+Q(NRAN+1)+E(NRAN))*X+Q(NRAN)*Q(NRAN+1)=0
   38 P=0.5D0*(Q(NRAN)+E(NRAN)+Q(NRAN+1))
      O=P*P-Q(NRAN)*Q(NRAN+1)
      T=DSQRT(DABS(O))
C
C        TEST FOR CONVERGENCE
      IF(S-TEPS)26,26,39
C
C        ARE THERE COMPLEX ROOTS
   39 IF(O)43,40,40
   40 IF(P)42,41,41
   41 T=-T
   42 P=P+T
      R=S
      GOTO 34
C
C        MODIFICATION FOR COMPLEX ROOTS
C        IS DISPLACEMENT SMALL ENOUGH
   43 IF(S-TDELT)44,35,35
C
C           INITIALIZATION
   44 O=Q(JBEG)+E(JBEG)-P
C
C           TEST FOR SMALL DIVISOR
      IF(DABS(O)-POL(JBEG))81,81,45
   45 T=(T/O)**2
      U=E(JBEG)*Q(JBEG+1)/(O*(1.0D0+T))
      V=O+U
C
C        THREEFOLD LOOP FOR COMPLEX DISPLACEMENT
      KOUNT=KOUNT+2
      DO 53 J=JBEG,NRAN
      O=Q(J+1)+E(J+1)-U-P
C
C           TEST FOR SMALL DIVISOR
      IF(DABS(V)-POL(J))46,46,49
   46 IF(J-NRAN)81,47,81
   47 EXPT=EXPT+P
      IF(ABS(SNGL(E(JEND)))-TOL)48,48,81
   48 P=0.5D0*(V+O-E(JEND))
      O=P*P-(V-U)*(O-U*T-O*W*(1.D0+T)/Q(JEND))
      T=DSQRT(DABS(O))
      GOTO 26
C
C           TEST FOR SMALL DIVISOR
   49 IF(DABS(O)-POL(J+1))46,46,50
   50 W=U*O/V
      T=T*(V/O)**2
      Q(J)=V+W-E(J-1)
      U=0.D0
      IF(J-NRAN)51,52,52
   51 U=Q(J+2)*E(J+1)/(O*(1.D0+T))
   52 V=O+U-W
C
C           TEST FOR SMALL DIVISOR
      IF(DABS(Q(J))-POL(J))81,81,53
   53 E(J)=W*V*(1.0D0+T)/Q(J)
      Q(NRAN+1)=V-E(NRAN)
   54 EXPT=EXPT+P
      TEPS=TEPS*1.1
      TDELT=TDELT*1.1
      IF(KOUNT-LIMIT)32,55,55
C
C        NO CONVERGENCE WITH FEASIBLE TOLERANCE
C           ERROR RETURN IN CASE OF UNSATISFACTORY CONVERGENCE
   55 IER=1
C
C        REARRANGE CALCULATED ROOTS
   56 IEND=NSAV-NRAN-1
      E(ISTA)=ESAV
      IF(IEND)59,59,57
   57 DO 58 I=1,IEND
      J=ISTA+I
      K=NRAN+1+I
      E(J)=E(K)
   58 Q(J)=Q(K)
   59 IR=ISTA+IEND
C
C        NORMAL RETURN
   60 IR=IR-1
      IF(IR)78,78,61
C
C        REARRANGE CALCULATED ROOTS
   61 DO 62 I=1,IR
      Q(I)=Q(I+1)
   62 E(I)=E(I+1)
C
C        CALCULATE COEFFICIENT VECTOR FROM ROOTS
      POL(IR+1)=1.D0
      IEND=IR-1
      JBEG=1
      DO 69 J=1,IR
      ISTA=IR+1-J
      O=0.D0
      P=Q(ISTA)
      T=E(ISTA)
      IF(T)65,63,65
C
C        MULTIPLY WITH LINEAR FACTOR
   63 DO 64 I=ISTA,IR
      POL(I)=O-P*POL(I+1)
   64 O=POL(I+1)
      GOTO 69
   65 GOTO(66,67),JBEG
   66 JBEG=2
      POL(ISTA)=0.D0
      GOTO 69
C
C        MULTIPLY WITH QUADRATIC FACTOR
   67 JBEG=1
      U=P*P+T*T
      P=P+P
      DO 68 I=ISTA,IEND
      POL(I)=O-P*POL(I+1)+U*POL(I+2)
   68 O=POL(I+1)
      POL(IR)=O-P
   69 CONTINUE
      IF(IER)78,70,78
C
C        COMPARISON OF COEFFICIENT VECTORS, IE. TEST OF ACCURACY
   70 P=0.D0
      DO 75 I=1,IR
      IF(C(I))72,71,72
   71 O=DABS(POL(I))
      GOTO 73
   72 O=DABS((POL(I)-C(I))/C(I))
   73 IF(P-O)74,75,75
   74 P=O
   75 CONTINUE
      IF(SNGL(P)-TOL)77,76,76
   76 IER=-1
   77 Q(IR+1)=P
      E(IR+1)=0.D0
   78 RETURN
C
C        ERROR RETURNS
C           ERROR RETURN FOR POLYNOMIALS OF DEGREE LESS THAN 1
   79 IER=2
      IR=0
      RETURN
C
C           ERROR RETURN IF THERE EXISTS NO S-FRACTION
   80 IER=4
      IR=ISTA
      GOTO 60
C
C           ERROR RETURN IN CASE OF INSTABLE QD-ALGORITHM
   81 IER=3
      GOTO 56
      END
